#pragma once
#include <Math/Matrix3x3.hpp>
#include <Math/Vector3.hpp>
namespace zzz{
template <std::size_t N,typename T=float>
struct SHCoeff
{
  T v[N*N];
  const int order_;
  SHCoeff():order_(N)
  {
  }
  const SHCoeff& operator=(float *data)
  {
    for (int i=0; i<N*N; i++) v[i]=data[i];
    return *this;
  }
  const SHCoeff& operator=(const SHCoeff<N,T>& coef)
  {
    for (int i=0; i<N*N; i++) v[i]=coef[i];
    return *this;
  }
  //NO +-*/ SINCE +=,-=,*=,/= IS RECOMMEND TO USE
  void operator+=(const SHCoeff<N,T>& coef)
  {
    for (int i=0; i<N*N; i++) v[i]+=coef.v[i];
  }
  void operator-=(const SHCoeff<N,T>& coef)
  {
    for (int i=0; i<N*N; i++) v[i]-=coef.v[i];
  }
  void operator*=(const T scaler)
  {
    for (int i=0; i<N*N; i++) v[i]*=scaler;
  }
  void operator/=(const T scaler)
  {
    for (int i=0; i<N*N; i++) v[i]/=scaler;
  }
  const SHCoeff operator*(float scale)
  {
    SHCoeff<N,T> ret;
    for (int i=0; i<N*N; i++) ret.v[i]=v[i]*scale;
    return ret;
  }
  inline T& operator [](int index)
  {
    return v[index];
  }
  inline const T& operator [](int index) const
  {
    return v[index];
  }
  float Dot(const SHCoeff<N,T> &coef)
  {
    float ret=0;
    for (int i=0; i<N*N; i++) ret+=v[i]*coef.v[i];
    return ret;
  }
  void Zero(T x=0)
  {
    for (int i=0; i<N*N; i++)
      v[i]=x;
  }
  void MultiMatRowMajor(const float *mat)
  {
    SHCoeff<N,T> tmp;
    for (int i=0; i<N*N; i++) for (int j=0; j<N*N; j++,mat++)
      tmp.v[i]+=*mat*v[j];
    *this=tmp;
  }
  void MultiMatColMajor(const float *mat)
  {
    SHCoeff<N,T> tmp;
    for (int i=0; i<N*N; i++) for (int j=0; j<N*N; j++,mat++)
      tmp.v[j]+=*mat*v[i];
    *this=tmp;
  }
  void MultiRotateMat(const float *mat)
  {
    SHCoeff<N,T> tmp;
    float *t=v;
    for (int i=0; i<N; i++)
      for (int j=0; j<=2*i; j++,mat+=N*N,t++)
        for (int k=i*i; k<(i+1)*(i+1); k++)
          *t+=mat[k]*v[k];
    *this=tmp;
  }
  void MultiRotateMatAtRightSize(const float *mat)
  {
    SHCoeff<N,T> tmp;
    float *t=v;
    for (int i=0; i<N; i++)
      for (int j=0; j<=2*i; j++,mat++,t++)
        for (int k=i*i; k<(i+1)*(i+1); k++)
          *t+=mat[k]*v[k];
    *this=tmp;
  }
  void TripleProduct(const SHCoeff<N,T> &coef)
  {
    printf("TripleProduct: Must be specialize to use!\n");
  }
  void TripleProductMatrix(float *ret)
  {
    printf("TripleProductMatrix: Must be specialize to use!\n");
  }
  friend ostream& operator<<(ostream &o,SHCoeff<N,T> &me)
  {
    me.SaveToFileA(o);
    return o;
  }
  friend istream& operator>>(istream &i,SHCoeff<N,T> &me)
  {
    me.LoadFromFileA(i);
    return i;
  }
  void SaveToFileA(ostream &fo)
  {
    for (int i=0; i<N*N; i++)
      fo<<v[i]<<' ';
    fo<<endl;
  }
  void LoadFromFileA(istream &fi)
  {
    for (int i=0; i<N*N; i++)
      fi>>v[i];
  }
  void LoadFromFileA(FILE *fp)
  {
    printf("LoadFromFileA(FILE *fp) must be specialized to use!\n");
  }
  void SaveToFileB(FILE *fp)
  {
    fwrite(v,sizeof(T),N*N,fp);
  }
  void LoadFromFileB(FILE *fp)
  {
    fread(v,sizeof(T),N*N,fp);
  }

//ROTATION
  void RotationMatrix_Ivanic(const Matrix<3,3,T> &rot)
  {
    int order=N;
    T ret[N*N*N*N];
    T R[10][10];
    T M[10][10*2-1][10*2-1];
    static T u[10][10*2-1][10*2-1],U[10][10*2-1][10*2-1],v[10][10*2-1][10*2-1],V[10][10*2-1][10*2-1],w[10][10*2-1][10*2-1],W[10][10*2-1][10*2-1];
    T P[10][3][10*2-1][10*2-1];
    memset(R,0,sizeof(R));
    memset(M,0,sizeof(M));
    memset(P,0,sizeof(P));
    static bool firsttime=true;
    if (firsttime)
    {
      memset(u,0,sizeof(u));
      memset(v,0,sizeof(v));
      memset(w,0,sizeof(w));
      firsttime=false;
    }
    memset(U,0,sizeof(U));
    memset(V,0,sizeof(V));
    memset(W,0,sizeof(W));

    R[0][0]=rot.v[4];  R[0][1]=-rot.v[7];  R[0][2]=rot.v[1];
    R[1][0]=-rot.v[5];  R[1][1]=rot.v[8];  R[1][2]=-rot.v[2];
    R[2][0]=rot.v[3];  R[2][1]=-rot.v[6];  R[2][2]=rot.v[0];

    M[0][0][0]=1.0f;
    M[1][0][0]=R[0][0];  M[1][0][1]=R[0][1];  M[1][0][2]=R[0][2];
    M[1][1][0]=R[1][0];  M[1][1][1]=R[1][1];  M[1][1][2]=R[1][2];
    M[1][2][0]=R[2][0];  M[1][2][1]=R[2][1];  M[1][2][2]=R[2][2];

    static int last_largest_order=2;

    if (order>last_largest_order)
    {
      for (int i=last_largest_order; i<order; i++)
        for (int s=-i; s<=i; s++)
          for (int t=-i; t<=i; t++)
          {
            if (abs(t)<i)
            {
              float d=(s==0?1.0f:0.0f);
              u[i][s+i][t+i]=sqrt((i+s)*(i-s)/(float)(i+t)/(i-t));
              v[i][s+i][t+i]=0.5f*sqrt((1+d)*(i+abs(s)-1)*(i+abs(s))/(float)(i+t)/(i-t))*(1.0f-2.0f*d);
              w[i][s+i][t+i]=-0.5f*sqrt((i-abs(s)-1)*(i-abs(s))/(float)(i+t)/(i-t))*(1.0f-d);
            }
            else
            {
              float d=(s==0?1.0f:0.0f);
              u[i][s+i][t+i]=sqrt((i+s)*(i-s)/(float)(2.0f*i)/(2.0f*i-1.0f));
              v[i][s+i][t+i]=0.5f*sqrt((1+d)*(i+abs(s)-1)*(i+abs(s))/(float)(2.0f*i)/(2.0f*i-1.0f))*(1.0f-2.0f*d);
              w[i][s+i][t+i]=-0.5f*sqrt((i-abs(s)-1)*(i-abs(s))/(float)(2.0f*i)/(2.0f*i-1.0f))*(1.0f-d);
            }
          }
          last_largest_order=order;
    }
    for (int i=2; i<order; i++)
    {
      for (int s=-i+1; s<=i-1; s++)
        for (int t=-i; t<=i; t++)
          for (int r=-1; r<=1; r++)
          {
            if (abs(t)<i)
              P[i][r+1][s+i][t+i]=R[r+1][1]*M[i-1][s+i-1][t+i-1];
            if (t==i)
              P[i][r+1][s+i][t+i]=R[r+1][2]*M[i-1][s+i-1][i*2-2]-R[r+1][0]*M[i-1][s+i-1][0];
            if (t==-i)
              P[i][r+1][s+i][t+i]=R[r+1][2]*M[i-1][s+i-1][0]+R[r+1][0]*M[i-1][s+i-1][i*2-2];
          }
          for (int t=-i; t<=i; t++)
          {
            for (int s=-i+1; s<=i-1; s++)
              U[i][s+i][t+i]=P[i][1][s+i][t+i];
            for (int s=-i; s<=i; s++)
            {
              if (s==0)
                V[i][s+i][t+i]=P[i][2][1+i][t+i]+P[i][0][-1+i][t+i];
              if (s>0)
              {
                float d=(s==1?1.0f:0.0f);
                V[i][s+i][t+i]=P[i][2][s-1+i][t+i]*sqrt(1.0f+d)-P[i][0][-s+1+i][t+i]*(1.0f-d);
              }
              if (s<0)
              {
                float d=(s==-1?1.0f:0.0f);
                V[i][s+i][t+i]=P[i][2][s+1+i][t+i]*(1.0f-d)+P[i][0][-s-1+i][t+i]*sqrt(1.0f+d);
              }
            }
            for (int s=-i+2; s<=i-2; s++)
            {
              if (s>0)
                W[i][s+i][t+i]=P[i][2][s+1+i][t+i]+P[i][0][-s-1+i][t+i];
              if (s<0)
                W[i][s+i][t+i]=P[i][2][s-1+i][t+i]-P[i][0][-s+1+i][t+i];
            }
          }
          for (int s=0; s<=i*2; s++)
            for (int t=0; t<=i*2; t++)
              M[i][s][t]=u[i][s][t]*U[i][s][t]+v[i][s][t]*V[i][s][t]+w[i][s][t]*W[i][s][t];
    }

    memset(ret,0,sizeof(float)*order*order*order*order);
    for (int i=0; i<order; i++)
      for (int j=0; j<=i*2; j++)
        for (int k=0; k<=i*2; k++)
          ret[(i*i+j)*order*order+i*i+k]=M[i][j][k];
    
    T newv[N*N];
    memset(newv,0,sizeof(T)*N*N);
    for (int i=0; i<N*N; i++)
      for (int j=0; j<N*N; j++)
        newv[i]+=this->v[j]*ret[j*N*N+i];
    memcpy(this->v,newv,sizeof(T)*N*N);

  }

  void Rotate_ZXZXZ(const float *rot)
  {
    printf("Rotate_ZXZXZ() must be specified to use!\n");
  }

  void MultRotateMat(const float *mat)        //time complexity 4/3*(order_)^3-(order_)/3
  {
    SHCoeff<N,T> ret;
    float *t=ret.v;
    for (int i=0; i<order_; i++)
      for (int j=0; j<=2*i; j++,mat+=order_*order_,t++)
        for (int k=i*i; k<(i+1)*(i+1); k++)
          *t+=mat[k]*v[k];
    *this=ret;
  }

  void MultRotateMatAtRightSide(const float *mat)
  {
    SHCoeff<N,T> ret;
    float *t=ret.v;
    for (int i=0; i<order_; i++)
      for (int j=0; j<=2*i; j++,mat++,t++)
        for (int k=i*i; k<(i+1)*(i+1); k++)
          *t+=mat[k*order_*order_]*v[k];
  }

};

template<>
void SHCoeff<4,float>::TripleProduct(const SHCoeff<4,float> &coef)
{
  SHCoeff<4,float> ret;
  float *c=ret.v;
  memset(c,0,sizeof(float)*16);
  const float *a=v,*b=coef.v;
  c[6]+=0.090109f*a[5]*b[5];
  c[8]+=-0.058340f*(a[3]*b[13]+a[13]*b[3]);
  c[0]+=0.281264f*a[9]*b[9];
  c[7]+=0.059470f*(a[12]*b[13]+a[13]*b[12]);
  c[5]+=0.059471f*(a[11]*b[12]+a[12]*b[11]);
  c[3]+=-0.058340f*(a[8]*b[13]+a[13]*b[8]);
  c[13]+=-0.058340f*(a[3]*b[8]+a[8]*b[3]);
  c[11]+=0.058369f*(a[1]*b[8]+a[8]*b[1]);
  c[1]+=0.058369f*(a[8]*b[11]+a[11]*b[8]);
  c[4]+=-0.058452f*(a[3]*b[11]+a[11]*b[3]);
  c[12]+=0.059470f*(a[7]*b[13]+a[13]*b[7]);
  c[6]+=0.090120f*a[7]*b[7];
  c[8]+=0.058369f*(a[1]*b[11]+a[11]*b[1]);
  c[0]+=0.281341f*a[8]*b[8];
  c[7]+=0.090120f*(a[6]*b[7]+a[7]*b[6]);
  c[5]+=0.090109f*(a[5]*b[6]+a[6]*b[5]);
  c[3]+=-0.058452f*(a[4]*b[11]+a[11]*b[4]);
  c[13]+=-0.058452f*(a[1]*b[4]+a[4]*b[1]);
  c[11]+=-0.058452f*(a[3]*b[4]+a[4]*b[3]);
  c[1]+=-0.058452f*(a[4]*b[13]+a[13]*b[4]);
  c[4]+=-0.058452f*(a[1]*b[13]+a[13]*b[1]);
  c[10]+=0.115150f*(a[5]*b[13]+a[13]*b[5]);
  c[6]+=0.126182f*a[11]*b[11];
  c[8]+=-0.094025f*(a[9]*b[11]+a[11]*b[9]);
  c[0]+=0.281361f*a[15]*b[15];
  c[7]+=0.115150f*(a[10]*b[11]+a[11]*b[10]);
  c[5]+=0.115150f*(a[10]*b[13]+a[13]*b[10]);
  c[3]+=-0.126223f*(a[3]*b[6]+a[6]*b[3]);
  c[13]+=0.059470f*(a[7]*b[12]+a[12]*b[7]);
  c[11]+=0.059471f*(a[5]*b[12]+a[12]*b[5]);
  c[1]+=-0.126236f*(a[1]*b[6]+a[6]*b[1]);
  c[4]+=-0.094033f*(a[9]*b[13]+a[13]*b[9]);
  c[14]+=0.115169f*(a[7]*b[13]+a[13]*b[7]);
  c[6]+=-0.126223f*a[3]*b[3];
  c[8]+=-0.094056f*(a[13]*b[15]+a[15]*b[13]);
  c[0]+=0.281660f*a[3]*b[3];
  c[7]+=0.115169f*(a[13]*b[14]+a[14]*b[13]);
  c[5]+=-0.115175f*(a[11]*b[14]+a[14]*b[11]);
  c[3]+=-0.142912f*(a[7]*b[12]+a[12]*b[7]);
  c[13]+=-0.094033f*(a[4]*b[9]+a[9]*b[4]);
  c[11]+=-0.094025f*(a[8]*b[9]+a[9]*b[8]);
  c[1]+=-0.143046f*(a[5]*b[12]+a[12]*b[5]);
  c[4]+=0.094033f*(a[11]*b[15]+a[15]*b[11]);
  c[2]+=0.184449f*(a[8]*b[14]+a[14]*b[8]);
  c[6]+=-0.126236f*a[1]*b[1];
  c[8]+=-0.145674f*a[11]*b[11];
  c[0]+=0.281727f*a[14]*b[14];
  c[7]+=-0.142912f*(a[3]*b[12]+a[12]*b[3]);
  c[5]+=-0.143046f*(a[1]*b[12]+a[12]*b[1]);
  c[3]+=0.184557f*(a[5]*b[10]+a[10]*b[5]);
  c[13]+=-0.094056f*(a[8]*b[15]+a[15]*b[8]);
  c[11]+=0.094033f*(a[4]*b[15]+a[15]*b[4]);
  c[1]+=0.184557f*(a[7]*b[10]+a[10]*b[7]);
  c[4]+=0.145561f*(a[11]*b[13]+a[13]*b[11]);
  c[12]+=0.059471f*(a[5]*b[11]+a[11]*b[5]);
  c[6]+=0.126282f*a[13]*b[13];
  c[8]+=0.145785f*a[13]*b[13];
  c[0]+=0.281767f*a[1]*b[1];
  c[7]+=0.148543f*(a[14]*b[15]+a[15]*b[14]);
  c[5]+=0.148600f*(a[9]*b[14]+a[14]*b[9]);
  c[9]+=-0.094025f*(a[8]*b[11]+a[11]*b[8]);
  c[15]+=0.094033f*(a[4]*b[11]+a[11]*b[4]);
  c[10]+=0.115150f*(a[7]*b[11]+a[11]*b[7]);
  c[11]+=0.115150f*(a[7]*b[10]+a[10]*b[7]);
  c[13]+=0.115150f*(a[5]*b[10]+a[10]*b[5]);
  c[14]+=-0.115175f*(a[5]*b[11]+a[11]*b[5]);
  c[6]+=0.168126f*a[12]*b[12];
  c[8]+=-0.156112f*a[5]*b[5];
  c[0]+=0.281884f*a[12]*b[12];
  c[7]+=0.148683f*(a[9]*b[10]+a[10]*b[9]);
  c[5]+=-0.148673f*(a[10]*b[15]+a[15]*b[10]);
  c[4]+=0.155967f*(a[5]*b[7]+a[7]*b[5]);
  c[2]+=0.184558f*(a[4]*b[10]+a[10]*b[4]);
  c[1]+=-0.184781f*(a[5]*b[14]+a[14]*b[5]);
  c[3]+=0.184955f*(a[7]*b[14]+a[14]*b[7]);
  c[9]+=-0.094033f*(a[4]*b[13]+a[13]*b[4]);
  c[15]+=-0.094056f*(a[8]*b[13]+a[13]*b[8]);
  c[6]+=0.180199f*a[6]*b[6];
  c[8]+=0.156254f*a[7]*b[7];
  c[0]+=0.281894f*a[10]*b[10];
  c[13]+=0.115169f*(a[7]*b[14]+a[14]*b[7]);
  c[11]+=-0.115175f*(a[5]*b[14]+a[14]*b[5]);
  c[12]+=-0.142912f*(a[3]*b[7]+a[7]*b[3]);
  c[14]+=0.148543f*(a[7]*b[15]+a[15]*b[7]);
  c[10]+=-0.148673f*(a[5]*b[15]+a[15]*b[5]);
  c[5]+=0.155967f*(a[4]*b[7]+a[7]*b[4]);
  c[7]+=0.155967f*(a[4]*b[5]+a[5]*b[4]);
  c[4]+=-0.180333f*(a[4]*b[6]+a[6]*b[4]);
  c[6]+=-0.180333f*a[4]*b[4];
  c[8]+=-0.180343f*(a[6]*b[8]+a[8]*b[6]);
  c[0]+=0.281897f*a[4]*b[4];
  c[3]+=0.202190f*(a[6]*b[13]+a[13]*b[6]);
  c[1]+=0.202243f*(a[6]*b[11]+a[11]*b[6]);
  c[2]+=0.218414f*(a[3]*b[7]+a[7]*b[3]);
  c[11]+=0.126182f*(a[6]*b[11]+a[11]*b[6]);
  c[13]+=0.126282f*(a[6]*b[13]+a[13]*b[6]);
  c[12]+=-0.143046f*(a[1]*b[5]+a[5]*b[1]);
  c[15]+=0.148543f*(a[7]*b[14]+a[14]*b[7]);
  c[14]+=0.148600f*(a[5]*b[9]+a[9]*b[5]);
  c[6]+=-0.180343f*a[8]*b[8];
  c[8]+=0.184449f*(a[2]*b[14]+a[14]*b[2]);
  c[0]+=0.281956f*a[13]*b[13];
  c[9]+=0.148600f*(a[5]*b[14]+a[14]*b[5]);
  c[10]+=0.148683f*(a[7]*b[9]+a[9]*b[7]);
  c[5]+=-0.156112f*(a[5]*b[8]+a[8]*b[5]);
  c[7]+=0.156254f*(a[7]*b[8]+a[8]*b[7]);
  c[4]+=0.184558f*(a[2]*b[10]+a[10]*b[2]);
  c[1]+=0.218374f*(a[3]*b[4]+a[4]*b[3]);
  c[3]+=0.218374f*(a[1]*b[4]+a[4]*b[1]);
  c[2]+=0.218461f*(a[1]*b[5]+a[5]*b[1]);
  c[6]+=0.202190f*(a[3]*b[13]+a[13]*b[3]);
  c[8]+=-0.187900f*(a[12]*b[14]+a[14]*b[12]);
  c[0]+=0.281966f*a[6]*b[6];
  c[11]+=0.145561f*(a[4]*b[13]+a[13]*b[4]);
  c[13]+=0.145561f*(a[4]*b[11]+a[11]*b[4]);
  c[15]+=-0.148673f*(a[5]*b[10]+a[10]*b[5]);
  c[9]+=0.148683f*(a[7]*b[10]+a[10]*b[7]);
  c[12]+=0.168126f*(a[6]*b[12]+a[12]*b[6]);
  c[14]+=0.184449f*(a[2]*b[8]+a[8]*b[2]);
  c[5]+=0.184557f*(a[3]*b[10]+a[10]*b[3]);
  c[10]+=0.184557f*(a[3]*b[5]+a[5]*b[3]);
  c[6]+=0.202243f*(a[1]*b[11]+a[11]*b[1]);
  c[7]+=0.184558f*(a[1]*b[10]+a[10]*b[1]);
  c[4]+=-0.188019f*(a[10]*b[12]+a[12]*b[10]);
  c[3]+=0.218414f*(a[2]*b[7]+a[7]*b[2]);
  c[1]+=0.218461f*(a[2]*b[5]+a[5]*b[2]);
  c[8]+=-0.218732f*a[1]*b[1];
  c[2]+=0.233572f*(a[5]*b[11]+a[11]*b[5]);
  c[0]+=0.281987f*a[7]*b[7];
  c[11]+=-0.145674f*(a[8]*b[11]+a[11]*b[8]);
  c[13]+=0.145785f*(a[8]*b[13]+a[13]*b[8]);
  c[10]+=0.184558f*(a[1]*b[7]+a[7]*b[1]);
  c[6]+=-0.210482f*a[9]*b[9];
  c[5]+=-0.184781f*(a[1]*b[14]+a[14]*b[1]);
  c[14]+=-0.184781f*(a[1]*b[5]+a[5]*b[1]);
  c[7]+=0.184955f*(a[3]*b[14]+a[14]*b[3]);
  c[12]+=-0.187900f*(a[8]*b[14]+a[14]*b[8]);
  c[9]+=-0.210482f*(a[6]*b[9]+a[9]*b[6]);
  c[15]+=-0.210522f*(a[6]*b[15]+a[15]*b[6]);
  c[4]+=0.218374f*(a[1]*b[3]+a[3]*b[1]);
  c[1]+=-0.218732f*(a[1]*b[8]+a[8]*b[1]);
  c[8]+=0.218821f*a[3]*b[3];
  c[3]+=0.218821f*(a[3]*b[8]+a[8]*b[3]);
  c[6]+=-0.210522f*a[15]*b[15];
  c[2]+=0.233583f*(a[7]*b[13]+a[13]*b[7]);
  c[0]+=0.282009f*a[5]*b[5];
  c[10]+=0.184558f*(a[2]*b[4]+a[4]*b[2]);
  c[14]+=0.184955f*(a[3]*b[7]+a[7]*b[3]);
  c[12]+=-0.188019f*(a[4]*b[10]+a[10]*b[4]);
  c[13]+=0.202190f*(a[3]*b[6]+a[6]*b[3]);
  c[11]+=0.202243f*(a[1]*b[6]+a[6]*b[1]);
  c[7]+=0.218414f*(a[2]*b[3]+a[3]*b[2]);
  c[5]+=0.218461f*(a[1]*b[2]+a[2]*b[1]);
  c[15]+=0.226034f*(a[3]*b[8]+a[8]*b[3]);
  c[3]+=0.226034f*(a[8]*b[15]+a[15]*b[8]);
  c[8]+=0.226034f*(a[3]*b[15]+a[15]*b[3]);
  c[1]+=0.226108f*(a[8]*b[9]+a[9]*b[8]);
  c[9]+=0.226108f*(a[1]*b[8]+a[8]*b[1]);
  c[4]+=-0.226136f*(a[1]*b[15]+a[15]*b[1]);
  c[6]+=0.247669f*(a[2]*b[12]+a[12]*b[2]);
  c[2]+=0.247669f*(a[6]*b[12]+a[12]*b[6]);
  c[0]+=0.282026f*a[11]*b[11];
  c[14]+=-0.187900f*(a[8]*b[12]+a[12]*b[8]);
  c[10]+=-0.188019f*(a[4]*b[12]+a[12]*b[4]);
  c[15]+=-0.226136f*(a[1]*b[4]+a[4]*b[1]);
  c[3]+=0.226185f*(a[4]*b[9]+a[9]*b[4]);
  c[8]+=0.226108f*(a[1]*b[9]+a[9]*b[1]);
  c[1]+=-0.226136f*(a[4]*b[15]+a[15]*b[4]);
  c[9]+=0.226185f*(a[3]*b[4]+a[4]*b[3]);
  c[4]+=0.226185f*(a[3]*b[9]+a[9]*b[3]);
  c[5]+=0.233572f*(a[2]*b[11]+a[11]*b[2]);
  c[11]+=0.233572f*(a[2]*b[5]+a[5]*b[2]);
  c[7]+=0.233583f*(a[2]*b[13]+a[13]*b[2]);
  c[13]+=0.233583f*(a[2]*b[7]+a[7]*b[2]);
  c[12]+=0.247669f*(a[2]*b[6]+a[6]*b[2]);
  c[2]+=0.252308f*(a[2]*b[6]+a[6]*b[2]);
  c[6]+=0.252308f*a[2]*b[2];
  c[0]+=0.282056f*a[0]*b[0];
  c[8]+=0.281341f*(a[0]*b[8]+a[8]*b[0]);
  c[9]+=0.281264f*(a[0]*b[9]+a[9]*b[0]);
  c[15]+=0.281361f*(a[0]*b[15]+a[15]*b[0]);
  c[3]+=0.281660f*(a[0]*b[3]+a[3]*b[0]);
  c[14]+=0.281727f*(a[0]*b[14]+a[14]*b[0]);
  c[1]+=0.281767f*(a[0]*b[1]+a[1]*b[0]);
  c[10]+=0.281894f*(a[0]*b[10]+a[10]*b[0]);
  c[12]+=0.281884f*(a[0]*b[12]+a[12]*b[0]);
  c[4]+=0.281897f*(a[0]*b[4]+a[4]*b[0]);
  c[13]+=0.281956f*(a[0]*b[13]+a[13]*b[0]);
  c[6]+=0.281966f*(a[0]*b[6]+a[6]*b[0]);
  c[7]+=0.281987f*(a[0]*b[7]+a[7]*b[0]);
  c[5]+=0.282009f*(a[0]*b[5]+a[5]*b[0]);
  c[11]+=0.282026f*(a[0]*b[11]+a[11]*b[0]);
  c[0]+=0.282066f*a[2]*b[2];
  c[2]+=0.282066f*(a[0]*b[2]+a[2]*b[0]);
  *this=ret;
}
template<>
void SHCoeff<4,float>::TripleProductMatrix(float *ret)
{
  memset(ret,0,sizeof(float)*256);
  float *c=ret;
  float *a=v;
  c[101]+=0.090109f*a[5];
  c[141]+=-0.058340f*a[3];
  c[131]+=-0.058340f*a[13];
  c[9]+=0.281264f*a[9];
  c[125]+=0.059470f*a[12];
  c[124]+=0.059470f*a[13];
  c[92]+=0.059471f*a[11];
  c[91]+=0.059471f*a[12];
  c[61]+=-0.058340f*a[8];
  c[56]+=-0.058340f*a[13];
  c[216]+=-0.058340f*a[3];
  c[211]+=-0.058340f*a[8];
  c[184]+=0.058369f*a[1];
  c[177]+=0.058369f*a[8];
  c[27]+=0.058369f*a[8];
  c[24]+=0.058369f*a[11];
  c[75]+=-0.058452f*a[3];
  c[67]+=-0.058452f*a[11];
  c[205]+=0.059470f*a[7];
  c[199]+=0.059470f*a[13];
  c[103]+=0.090120f*a[7];
  c[139]+=0.058369f*a[1];
  c[129]+=0.058369f*a[11];
  c[8]+=0.281341f*a[8];
  c[119]+=0.090120f*a[6];
  c[118]+=0.090120f*a[7];
  c[86]+=0.090109f*a[5];
  c[85]+=0.090109f*a[6];
  c[59]+=-0.058452f*a[4];
  c[52]+=-0.058452f*a[11];
  c[212]+=-0.058452f*a[1];
  c[209]+=-0.058452f*a[4];
  c[180]+=-0.058452f*a[3];
  c[179]+=-0.058452f*a[4];
  c[29]+=-0.058452f*a[4];
  c[20]+=-0.058452f*a[13];
  c[77]+=-0.058452f*a[1];
  c[65]+=-0.058452f*a[13];
  c[173]+=0.115150f*a[5];
  c[165]+=0.115150f*a[13];
  c[107]+=0.126182f*a[11];
  c[139]+=-0.094025f*a[9];
  c[137]+=-0.094025f*a[11];
  c[15]+=0.281361f*a[15];
  c[123]+=0.115150f*a[10];
  c[122]+=0.115150f*a[11];
  c[93]+=0.115150f*a[10];
  c[90]+=0.115150f*a[13];
  c[54]+=-0.126223f*a[3];
  c[51]+=-0.126223f*a[6];
  c[220]+=0.059470f*a[7];
  c[215]+=0.059470f*a[12];
  c[188]+=0.059471f*a[5];
  c[181]+=0.059471f*a[12];
  c[22]+=-0.126236f*a[1];
  c[17]+=-0.126236f*a[6];
  c[77]+=-0.094033f*a[9];
  c[73]+=-0.094033f*a[13];
  c[237]+=0.115169f*a[7];
  c[231]+=0.115169f*a[13];
  c[99]+=-0.126223f*a[3];
  c[143]+=-0.094056f*a[13];
  c[141]+=-0.094056f*a[15];
  c[3]+=0.281660f*a[3];
  c[126]+=0.115169f*a[13];
  c[125]+=0.115169f*a[14];
  c[94]+=-0.115175f*a[11];
  c[91]+=-0.115175f*a[14];
  c[60]+=-0.142912f*a[7];
  c[55]+=-0.142912f*a[12];
  c[217]+=-0.094033f*a[4];
  c[212]+=-0.094033f*a[9];
  c[185]+=-0.094025f*a[8];
  c[184]+=-0.094025f*a[9];
  c[28]+=-0.143046f*a[5];
  c[21]+=-0.143046f*a[12];
  c[79]+=0.094033f*a[11];
  c[75]+=0.094033f*a[15];
  c[46]+=0.184449f*a[8];
  c[40]+=0.184449f*a[14];
  c[97]+=-0.126236f*a[1];
  c[139]+=-0.145674f*a[11];
  c[14]+=0.281727f*a[14];
  c[124]+=-0.142912f*a[3];
  c[115]+=-0.142912f*a[12];
  c[92]+=-0.143046f*a[1];
  c[81]+=-0.143046f*a[12];
  c[58]+=0.184557f*a[5];
  c[53]+=0.184557f*a[10];
  c[223]+=-0.094056f*a[8];
  c[216]+=-0.094056f*a[15];
  c[191]+=0.094033f*a[4];
  c[180]+=0.094033f*a[15];
  c[26]+=0.184557f*a[7];
  c[23]+=0.184557f*a[10];
  c[77]+=0.145561f*a[11];
  c[75]+=0.145561f*a[13];
  c[203]+=0.059471f*a[5];
  c[197]+=0.059471f*a[11];
  c[109]+=0.126282f*a[13];
  c[141]+=0.145785f*a[13];
  c[1]+=0.281767f*a[1];
  c[127]+=0.148543f*a[14];
  c[126]+=0.148543f*a[15];
  c[94]+=0.148600f*a[9];
  c[89]+=0.148600f*a[14];
  c[155]+=-0.094025f*a[8];
  c[152]+=-0.094025f*a[11];
  c[251]+=0.094033f*a[4];
  c[244]+=0.094033f*a[11];
  c[171]+=0.115150f*a[7];
  c[167]+=0.115150f*a[11];
  c[186]+=0.115150f*a[7];
  c[183]+=0.115150f*a[10];
  c[218]+=0.115150f*a[5];
  c[213]+=0.115150f*a[10];
  c[235]+=-0.115175f*a[5];
  c[229]+=-0.115175f*a[11];
  c[108]+=0.168126f*a[12];
  c[133]+=-0.156112f*a[5];
  c[12]+=0.281884f*a[12];
  c[122]+=0.148683f*a[9];
  c[121]+=0.148683f*a[10];
  c[95]+=-0.148673f*a[10];
  c[90]+=-0.148673f*a[15];
  c[71]+=0.155967f*a[5];
  c[69]+=0.155967f*a[7];
  c[42]+=0.184558f*a[4];
  c[36]+=0.184558f*a[10];
  c[30]+=-0.184781f*a[5];
  c[21]+=-0.184781f*a[14];
  c[62]+=0.184955f*a[7];
  c[55]+=0.184955f*a[14];
  c[157]+=-0.094033f*a[4];
  c[148]+=-0.094033f*a[13];
  c[253]+=-0.094056f*a[8];
  c[248]+=-0.094056f*a[13];
  c[102]+=0.180199f*a[6];
  c[135]+=0.156254f*a[7];
  c[10]+=0.281894f*a[10];
  c[222]+=0.115169f*a[7];
  c[215]+=0.115169f*a[14];
  c[190]+=-0.115175f*a[5];
  c[181]+=-0.115175f*a[14];
  c[199]+=-0.142912f*a[3];
  c[195]+=-0.142912f*a[7];
  c[239]+=0.148543f*a[7];
  c[231]+=0.148543f*a[15];
  c[175]+=-0.148673f*a[5];
  c[165]+=-0.148673f*a[15];
  c[87]+=0.155967f*a[4];
  c[84]+=0.155967f*a[7];
  c[117]+=0.155967f*a[4];
  c[116]+=0.155967f*a[5];
  c[70]+=-0.180333f*a[4];
  c[68]+=-0.180333f*a[6];
  c[100]+=-0.180333f*a[4];
  c[136]+=-0.180343f*a[6];
  c[134]+=-0.180343f*a[8];
  c[4]+=0.281897f*a[4];
  c[61]+=0.202190f*a[6];
  c[54]+=0.202190f*a[13];
  c[27]+=0.202243f*a[6];
  c[22]+=0.202243f*a[11];
  c[39]+=0.218414f*a[3];
  c[35]+=0.218414f*a[7];
  c[187]+=0.126182f*a[6];
  c[182]+=0.126182f*a[11];
  c[221]+=0.126282f*a[6];
  c[214]+=0.126282f*a[13];
  c[197]+=-0.143046f*a[1];
  c[193]+=-0.143046f*a[5];
  c[254]+=0.148543f*a[7];
  c[247]+=0.148543f*a[14];
  c[233]+=0.148600f*a[5];
  c[229]+=0.148600f*a[9];
  c[104]+=-0.180343f*a[8];
  c[142]+=0.184449f*a[2];
  c[130]+=0.184449f*a[14];
  c[13]+=0.281956f*a[13];
  c[158]+=0.148600f*a[5];
  c[149]+=0.148600f*a[14];
  c[169]+=0.148683f*a[7];
  c[167]+=0.148683f*a[9];
  c[88]+=-0.156112f*a[5];
  c[85]+=-0.156112f*a[8];
  c[120]+=0.156254f*a[7];
  c[119]+=0.156254f*a[8];
  c[74]+=0.184558f*a[2];
  c[66]+=0.184558f*a[10];
  c[20]+=0.218374f*a[3];
  c[19]+=0.218374f*a[4];
  c[52]+=0.218374f*a[1];
  c[49]+=0.218374f*a[4];
  c[37]+=0.218461f*a[1];
  c[33]+=0.218461f*a[5];
  c[109]+=0.202190f*a[3];
  c[99]+=0.202190f*a[13];
  c[142]+=-0.187900f*a[12];
  c[140]+=-0.187900f*a[14];
  c[6]+=0.281966f*a[6];
  c[189]+=0.145561f*a[4];
  c[180]+=0.145561f*a[13];
  c[219]+=0.145561f*a[4];
  c[212]+=0.145561f*a[11];
  c[250]+=-0.148673f*a[5];
  c[245]+=-0.148673f*a[10];
  c[154]+=0.148683f*a[7];
  c[151]+=0.148683f*a[10];
  c[204]+=0.168126f*a[6];
  c[198]+=0.168126f*a[12];
  c[232]+=0.184449f*a[2];
  c[226]+=0.184449f*a[8];
  c[90]+=0.184557f*a[3];
  c[83]+=0.184557f*a[10];
  c[165]+=0.184557f*a[3];
  c[163]+=0.184557f*a[5];
  c[107]+=0.202243f*a[1];
  c[97]+=0.202243f*a[11];
  c[122]+=0.184558f*a[1];
  c[113]+=0.184558f*a[10];
  c[76]+=-0.188019f*a[10];
  c[74]+=-0.188019f*a[12];
  c[55]+=0.218414f*a[2];
  c[50]+=0.218414f*a[7];
  c[21]+=0.218461f*a[2];
  c[18]+=0.218461f*a[5];
  c[129]+=-0.218732f*a[1];
  c[43]+=0.233572f*a[5];
  c[37]+=0.233572f*a[11];
  c[7]+=0.281987f*a[7];
  c[187]+=-0.145674f*a[8];
  c[184]+=-0.145674f*a[11];
  c[221]+=0.145785f*a[8];
  c[216]+=0.145785f*a[13];
  c[167]+=0.184558f*a[1];
  c[161]+=0.184558f*a[7];
  c[105]+=-0.210482f*a[9];
  c[94]+=-0.184781f*a[1];
  c[81]+=-0.184781f*a[14];
  c[229]+=-0.184781f*a[1];
  c[225]+=-0.184781f*a[5];
  c[126]+=0.184955f*a[3];
  c[115]+=0.184955f*a[14];
  c[206]+=-0.187900f*a[8];
  c[200]+=-0.187900f*a[14];
  c[153]+=-0.210482f*a[6];
  c[150]+=-0.210482f*a[9];
  c[255]+=-0.210522f*a[6];
  c[246]+=-0.210522f*a[15];
  c[67]+=0.218374f*a[1];
  c[65]+=0.218374f*a[3];
  c[24]+=-0.218732f*a[1];
  c[17]+=-0.218732f*a[8];
  c[131]+=0.218821f*a[3];
  c[56]+=0.218821f*a[3];
  c[51]+=0.218821f*a[8];
  c[111]+=-0.210522f*a[15];
  c[45]+=0.233583f*a[7];
  c[39]+=0.233583f*a[13];
  c[5]+=0.282009f*a[5];
  c[164]+=0.184558f*a[2];
  c[162]+=0.184558f*a[4];
  c[231]+=0.184955f*a[3];
  c[227]+=0.184955f*a[7];
  c[202]+=-0.188019f*a[4];
  c[196]+=-0.188019f*a[10];
  c[214]+=0.202190f*a[3];
  c[211]+=0.202190f*a[6];
  c[182]+=0.202243f*a[1];
  c[177]+=0.202243f*a[6];
  c[115]+=0.218414f*a[2];
  c[114]+=0.218414f*a[3];
  c[82]+=0.218461f*a[1];
  c[81]+=0.218461f*a[2];
  c[248]+=0.226034f*a[3];
  c[243]+=0.226034f*a[8];
  c[63]+=0.226034f*a[8];
  c[56]+=0.226034f*a[15];
  c[143]+=0.226034f*a[3];
  c[131]+=0.226034f*a[15];
  c[25]+=0.226108f*a[8];
  c[24]+=0.226108f*a[9];
  c[152]+=0.226108f*a[1];
  c[145]+=0.226108f*a[8];
  c[79]+=-0.226136f*a[1];
  c[65]+=-0.226136f*a[15];
  c[108]+=0.247669f*a[2];
  c[98]+=0.247669f*a[12];
  c[44]+=0.247669f*a[6];
  c[38]+=0.247669f*a[12];
  c[11]+=0.282026f*a[11];
  c[236]+=-0.187900f*a[8];
  c[232]+=-0.187900f*a[12];
  c[172]+=-0.188019f*a[4];
  c[164]+=-0.188019f*a[12];
  c[244]+=-0.226136f*a[1];
  c[241]+=-0.226136f*a[4];
  c[57]+=0.226185f*a[4];
  c[52]+=0.226185f*a[9];
  c[137]+=0.226108f*a[1];
  c[129]+=0.226108f*a[9];
  c[31]+=-0.226136f*a[4];
  c[20]+=-0.226136f*a[15];
  c[148]+=0.226185f*a[3];
  c[147]+=0.226185f*a[4];
  c[73]+=0.226185f*a[3];
  c[67]+=0.226185f*a[9];
  c[91]+=0.233572f*a[2];
  c[82]+=0.233572f*a[11];
  c[181]+=0.233572f*a[2];
  c[178]+=0.233572f*a[5];
  c[125]+=0.233583f*a[2];
  c[114]+=0.233583f*a[13];
  c[215]+=0.233583f*a[2];
  c[210]+=0.233583f*a[7];
  c[198]+=0.247669f*a[2];
  c[194]+=0.247669f*a[6];
  c[38]+=0.252308f*a[2];
  c[34]+=0.252308f*a[6];
  c[98]+=0.252308f*a[2];
  c[0]+=0.282056f*a[0];
  c[136]+=0.281341f*a[0];
  c[128]+=0.281341f*a[8];
  c[153]+=0.281264f*a[0];
  c[144]+=0.281264f*a[9];
  c[255]+=0.281361f*a[0];
  c[240]+=0.281361f*a[15];
  c[51]+=0.281660f*a[0];
  c[48]+=0.281660f*a[3];
  c[238]+=0.281727f*a[0];
  c[224]+=0.281727f*a[14];
  c[17]+=0.281767f*a[0];
  c[16]+=0.281767f*a[1];
  c[170]+=0.281894f*a[0];
  c[160]+=0.281894f*a[10];
  c[204]+=0.281884f*a[0];
  c[192]+=0.281884f*a[12];
  c[68]+=0.281897f*a[0];
  c[64]+=0.281897f*a[4];
  c[221]+=0.281956f*a[0];
  c[208]+=0.281956f*a[13];
  c[102]+=0.281966f*a[0];
  c[96]+=0.281966f*a[6];
  c[119]+=0.281987f*a[0];
  c[112]+=0.281987f*a[7];
  c[85]+=0.282009f*a[0];
  c[80]+=0.282009f*a[5];
  c[187]+=0.282026f*a[0];
  c[176]+=0.282026f*a[11];
  c[2]+=0.282066f*a[2];
  c[34]+=0.282066f*a[0];
  c[32]+=0.282066f*a[2];
}
template<>
void SHCoeff<4,float>::LoadFromFileA(FILE *fp)
{
  for (int i=0; i<4*4; i++)
    fscanf(fp,"%f",v+i);
}
template<>
void SHCoeff<4,Vector<3,float> >::LoadFromFileA(FILE *fp)
{
  for (int i=0; i<4*4; i++)
    fscanf(fp,"%f %f %f",&(v[i][0]),&(v[i][1]),&(v[i][2]));
}
template<>
void SHCoeff<4,float>::Rotate_ZXZXZ(const float *rot)
{
  SHCoeff<4,float> tmp(*this);

  float cosine[4];
  float sine[4];
#define CSHCOEF4F_ROTATE4_ZXZXZ_ROTATEZ(tar,src,v_cos,v_sin)\
  {\
  cosine[1]=v_cos; \
  sine[1]=v_sin; \
  cosine[2]=2.0f*cosine[1]*cosine[1]-1.0f; \
  sine[2]=2.0f*cosine[1]*sine[1]; \
  cosine[3]=cosine[2]*cosine[1]+sine[2]*sine[1]; \
  sine[3]=cosine[2]*sine[1]+sine[2]*cosine[1]; \
  tar[0]=src[0]; \
  tar[1]=cosine[1]*src[1]+sine[1]*src[3]; \
  tar[2]=src[2]; \
  tar[3]=cosine[1]*src[3]-sine[1]*src[1]; \
  tar[4]=cosine[2]*src[4]+sine[2]*src[8]; \
  tar[5]=cosine[1]*src[5]+sine[1]*src[7]; \
  tar[6]=src[6]; \
  tar[7]=cosine[1]*src[7]-sine[1]*src[5]; \
  tar[8]=cosine[2]*src[8]-sine[2]*src[4]; \
  tar[9]=cosine[3]*src[9]+sine[3]*src[15]; \
  tar[10]=cosine[2]*src[10]+sine[2]*src[14]; \
  tar[11]=cosine[1]*src[11]+sine[1]*src[13]; \
  tar[12]=src[12]; \
  tar[13]=cosine[1]*src[13]-sine[1]*src[11]; \
  tar[14]=cosine[2]*src[14]-sine[2]*src[10]; \
  tar[15]=cosine[3]*src[15]-sine[3]*src[9]; \
  }
  //need for SSE2
  /*_declspec(align(16))*/ float buf0[16];
  /*_declspec(align(16))*/ float buf1[16];

  if (1.0-rot[8]*rot[8]<1e-4 && rot[8]>0.0f)
  {
    CSHCOEF4F_ROTATE4_ZXZXZ_ROTATEZ(buf0,tmp.v,rot[0],-rot[3]);
    memcpy(v,buf0,sizeof(float)*16);
    return;
  }

  float cosa,sina,cosb,sinb,cosc,sinc;
  cosb=rot[8];
  sinb=sqrt(1-cosb*cosb);
  cosa=rot[2]/sinb;
  sina=rot[5]/sinb;
  cosc=-rot[6]/sinb;
  sinc=rot[7]/sinb;
  if (1.0-rot[8]*rot[8]<1e-4 && rot[8]<0.0f)
  {
    cosa=1.0f;
    sina=0.0f;
    cosc=-rot[0];
    sinc=rot[3];
  }

  CSHCOEF4F_ROTATE4_ZXZXZ_ROTATEZ(buf0,tmp.v,cosc,sinc);

  buf1[0]=buf0[0];
  buf1[1]=buf0[2];
  buf1[2]=-buf0[1];
  buf1[3]=buf0[3];
  buf1[4]=buf0[7];
  buf1[5]=-buf0[5];
  buf1[6]=-0.5f*buf0[6]-0.86603f*buf0[8];
  buf1[7]=-buf0[4];
  buf1[8]=-0.86603f*buf0[6]+0.5f*buf0[8];
  buf1[9]=-0.79057f*buf0[12]+0.61237f*buf0[14];
  buf1[10]=-buf0[10];
  buf1[11]=buf0[12]*-0.61237f+buf0[14]*-0.79057f;
  buf1[12]=buf0[9]*0.79057f+buf0[11]*0.61237f;
  buf1[13]=buf0[13]*-0.25f+buf0[15]*-0.96825f;
  buf1[14]=buf0[9]*-0.61237f+buf0[11]*0.79057f;
  buf1[15]=buf0[13]*-0.96825f+buf0[15]*0.25f;


  CSHCOEF4F_ROTATE4_ZXZXZ_ROTATEZ(buf0,buf1,cosb,sinb);

  buf1[0]=buf0[0];
  buf1[1]=-buf0[2];
  buf1[2]=buf0[1];
  buf1[3]=buf0[3];
  buf1[4]=-buf0[7];
  buf1[5]=-buf0[5];
  buf1[6]=buf0[6]*-0.5f+buf0[8]*-0.86603f;
  buf1[7]=buf0[4];
  buf1[8]=buf0[6]*-0.86603f+buf0[8]*0.5f;
  buf1[9]=buf0[12]*0.79057f+buf0[14]*-0.61237f;
  buf1[10]=-buf0[10];
  buf1[11]=buf0[12]*0.61237f+buf0[14]*0.79057f;
  buf1[12]=buf0[9]*-0.79057f+buf0[11]*-0.61237f;
  buf1[13]=buf0[13]*-0.25f+buf0[15]*-0.96825f;
  buf1[14]=buf0[9]*0.61237f+buf0[11]*-0.79057f;
  buf1[15]=buf0[13]*-0.96825f+buf0[15]*0.25f;

  CSHCOEF4F_ROTATE4_ZXZXZ_ROTATEZ(buf0,buf1,cosa,sina);

  memcpy(v,buf0,sizeof(float)*16);
#undef CSHCOEF4F_ROTATE4_ZXZXZ_ROTATEZ
}

typedef SHCoeff<4,float> SHCoeff4f;
typedef SHCoeff<4,Vector3f> SHCoeff4v3f;

}
